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[i] Electrostatic solitary structures have been observed in the Earth’s magnetosheath by 
the Cluster spacecraft. Recent theoretical work has suggested that these solitary structures 
are modeled by electron acoustic solitary waves existing in a four-component plasma 
system consisting of core electrons, two counter-streaming electron beams, and one species 
of background ions. In this paper, the excitation of electron acoustic waves and the 
formation of solitary structures are studied by means of a one-dimensional electrostatic 
Vlasov simulation. The present result first shows that either electron acoustic solitary 
waves with negative potential or electron phase-space holes with positive potential are 
excited in four-component plasma systems. However, these electrostatic solitary structures 
have longer duration times and higher wave amplitudes than the solitary structures 
observed in the magnetosheath. The result indicates that a high-speed and small free energy 
source may be needed as a fifth component. An additional simulation of a five-component 
plasma consisting of a stable four-component plasma and a weak electron beam shows 
the generation of small and fast electron phase-space holes by the bump-on-tail instability. 

The physical properties of the small and fast electron phase-space holes are very similar to 
those obtained by the previous theoretical analysis. The amplitude and duration time of 
solitary structures in the simulation are also in agreement with the Cluster observation. 

Citation: Umeda, T., M. Ashour-Abdalla, J. S. Pickett, and M. L. Goldstein (2012), Vlasov simulation of electrostatic solitary 
structures in multi-component plasmas, J. Geophys. Res., 117, A05223, doi: 10. 1029/201 1 JAO 17181. 


1. Introduction 

[ 2 ] The existence of Electrostatic Solitary Waves (ESWs) 
was first finnly established through observations made by the 
GEOTAIL spacecraft in the Earth’s magnetotail [Matsumoto 
etal., 1994]. These solitary waves have waveforms of bipolar 
electric field pulses longitudinal to the ambient magnetic 
field [Kojima et al., 1997], which suggests that the solitary 
waves are electrostatic solitary structures. The electrostatic 
solitary structures are also observed in other regions of the 
Earth’s magnetosphere: at the bow shock [Bale et al., 1998], 
in the auroral region [Mozer et al., 1997; Ergun et al., 1998], 
in the polar cap boundary layer [Franz et al., 1998], in the 
magnetotail magnetic reconnection region [Cattell et al., 
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2005], at the magnetopause [Cattell et al., 2002], and in the 
magnetosheath [Pickett et al., 2003, 2005]. 

[ 3 ] There are two basic physical models proposed for the 
generation of electrostatic solitary structures. One of these 
interprets the solitary structures as electron phase-space holes 
and the other as electron acoustic solitary waves. Previous 
computer simulation studies demonstrated that electron phase- 
space holes are generated during the nonlinear evolution of 
weak-electron-beam (bump-on-tail) instabilities [Omura et al. , 
1996; Miyake et al., 1998; Singh et al., 2000; Umeda et al., 
2002, 2004] and strong-electron-beam (two-stream) instabil- 
ities [Goldman et al., 1999; Miyake et al., 2000; Oppenheim 
et al., 2001; Umeda et al., 2006; Umeda, 2008a]. The elec- 
tron phase-space holes correspond to positive potentials that 
trap electrons. They are modeled by using one-dimensional 
equilibrium solutions to the time-independent Vlasov-Poisson 
equations, [Krasovsky et al., 1997, 2003; Chen and Parks, 
2002; Chen et al, 2004] which are called the Bemstein- 
Greene-Kruskal (BGK) equilibrium [Bernstein et al., 1957]. 

[ 4 ] Recently, there has been an attempt to explain the elec- 
trostatic solitary structures as density structures (enhancements 
or decreases) which result from ion or electron acoustic 
instabilities [e.g., Ghosh et al., 2008; Lakhina et al., 2008, 
2009, and references therein]. In the theoretical models of 
electron acoustic solitary waves, multicomponent plasma 
systems consisting of cold and hot electrons and one species 
of ions are assumed, and multifluid equations for each 
species are solved in a moving stationary frame. Lakhina 
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Table 1 . Parameters for Different Simulation Runs 


Run 

k' pee 

V dee 

v tce 

^ppe 

V dpe 

V tpe 

^pae 

Vdae 

v tae 

Upi 

v di 

Vti 

(dpbe 

V dbe 

V, he 

0 

0.7484 

0.0042 

1.0 

0.5023 

1.3152 

0.6685 

0.4331 

-1.3804 

0.5435 

0.0025 

0.0 

0.0344 

— 

— 

— 

1 

0.7484 

0.0042 

1.0 

0.5023 

1.3152 

0.4727 

0.4331 

-1.3804 

0.3843 

0.0025 

0.0 

0.0344 

- 

- 

- 

2 

0.7484 

0.0042 

1.0 

0.5023 

1.3152 

0.3342 

0.4331 

-1.3804 

0.2717 

0.0025 

0.0 

0.0344 

- 

- 

- 

3 

0.7484 

0.0042 

1.0 

0.5023 

1.3152 

0.6685 

0.4331 

-1.3804 

0.5435 

0.0025 

0.0 

0.0344 

0.1 

-2.2 

0.1 


et al. [2009] obtained solutions of positive solitary potentials 
by assuming four-component plasmas based on the Cluster 
observations in the magnetosheath [Pickett et al., 2005]. The 
electron parameters, in particular, were obtained from 
detailed analysis of electron distribution functions obtained 
by the PEACE (Plasma Electron And Current Experiment) 
instrument on-board one of the four Cluster spacecraft (SC4) 
[. Johnstone et al., 1997]. These distribution functions are 
obtained, at best, every 2 s, which is not atypical for most 
space-based particle measurements. Considering that typical 
instabilities, such as the acoustic instability that can generate 
solitary waves of short time duration usually develop over 
only hundreds to thousands of plasma periods (order of 10s 
of ms or less for this specific case), it is not clear that those 
distribution functions provide the pre-unstable parameters 
necessary for carrying out an investigation. Certainly, these 
distribution functions will represent a composite of pre- and 
post-instability characteristics, assuming that the solitary 
waves were, in fact, generated during the time of the particle 
measurement. Given that solitary waves observed in the 
magnetosheath can propagate over distances at least as large 
as a km [Pickett et al, 2010], it is not necessarily given that 
solitary waves observed in space are generated at the time 
and place of the measurement. However, it seems quite 
plausible to assume that they may be generated locally when 
observations show solitary waves of similar or sequential 
increasing/decreasing amplitudes spaced at discrete time 
intervals, such as the case well chosen by Lakhina et al. 
[2009, Figure 1]. 

[5] In spite of some of the measurement limitations and 
the lack of certainty of whether the solitary waves are gen- 
erated locally, the theoretical studies of the acoustic model 
were extremely useful in showing that there exist steady 
state solutions of electrostatic solitary structures in four- 
component plasmas. We now take these studies one step 
further by examining the excitation of electrostatic waves in 
four-component plasmas and their nonlinear evolution by 
means of a self-consistent solution of the Vlasov-Poisson 
equations. We will use observed input parameters initially 
and then vary some of those parameters slightly in order to 
explore the possibility that we can infer the state of the 
unperturbed plasma and the resultant excitation that pro- 
duces the observed solitary waves. 

[6] This paper is structured as follows. Section 2 provides 
the model and the parameters used in the computer simula- 
tions. Section 3 presents the simulation results. Section 4 
contains a summary of the computer simulation results 
together with a comparison of those results with spacecraft 
observations. 

2. Simulation Model and Parameters 

[7] In the present study, we solve spatiotemporal devel- 
opment of the complete position-velocity distribution 


function self-consistently based on the Vlasov equation 
together with the Poisson equation. The present numerical 
simulation code is a general lxl v (one position dimension 
and one velocity dimension) Vlasov code with the ambient 
magnetic field taken along the simulation domain. We 
ignore electromagnetic perturbations so that the effects of 
the ambient magnetic field and electromagnetic waves are 
completely absent in the simulation domain. Then, the dis- 
tribution function is reduced to the two-dimensional posi- 
tion-velocity phase space fix, v x ). The Vlasov simulation has 
significant benefit because numerical noise is suppressed 
substantially in comparison with one-dimensional electro- 
static particle-in-cell simulations. The Vlasov code adopts 
the time-splitting algorithm [Cheng and Knorr, 1976] and 
a non-oscillatory, positive and conservative interpolation 
scheme for stable time-integration of phase-space distribu- 
tion functions [Umeda, 2008b; Umeda et al., 2012], 

[8] Velocity distribution functions of the magnetosheath 
plasma are modeled first with four Maxwellian components 
[Lakhina et al., 2009]. 


f(Vx) 


V2nVts 


exp 


(V.r - Vdsf 

iy2 
/LV ts 


(1) 


These components are core electrons (the subscript s = ce), 
an electron beam propagating parallel to the magnetic field 
{s = pe), an electron beam propagating antiparallel to the 
magnetic field (s = ae), and background ions (5 = i). Simu- 
lation parameters for the present study are based on those 
listed in the paper by Lakhina et al. [2009], which were 
observed with the PEACE instrument [Johnstone et al., 
1997]. The normalized simulation parameters are listed in 
Table 1. Run 0 corresponds to “event 1” in the paper by 
Lakhina et al. [2009, Table 1]. Here, the plasma frequency 
(10 p S = e 2 njc( ) m s ) is normalized to the total electron plasma 
frequency (co pce + ca ppe + uj 2 pae = io 2 pe = 1), and the velocity is 
normalized to the thermal velocity of core electrons V tec . The 
ions are assumed to have a realistic mass ratio mjm e = 1600 
with thermal energy of 100 eV. 

[ 9 ] Although the initial velocity distribution function for 
Run 0 (i.e., the observed velocity distribution function) has a 
positive slope as shown by the dashed line in Figure 1 
(bottom), the linear dispersion relation does not give unsta- 
ble electron acoustic wave modes but rather gives only a 
stable Langmuir wave mode. There is also no wave growth 
in the Vlasov simulation for Run 0 (not shown). However, 
this is not surprising because a spacecraft sometimes 
observes distribution functions after saturation of micro- 
scopic instabilities and timescales for the saturation of elec- 
tron-scale microscopic instabilities are sometimes shorter 
than the time resolution of particle instruments. In order to 
make the initial distribution functions unstable, we need to 
decrease the electron temperature. 
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Figure 1 . The electron velocity distribution functions at the 
initial and final (io pe t = 4000) states for Runs 1-3. The elec- 
tron velocity distribution function for Run 0 is also plotted 
with dashed lines for reference. 


cold electrons at rest [e.g., Ashour-Abdalla and Okuda, 
1986]. In the present study, on the other hand, there are 
two counter- streaming cold electron beams and one hot 
electron component at rest, which is a nonstandard situation 
for electron acoustic waves. It is obvious that unstable 
modes are neither electron beam modes nor ion acoustic 
modes because the phase velocities are much slower than the 
drift velocities of the two electron beams and the core elec- 
tron temperature (45.8 eV) is smaller than the ion tempera- 
ture. The phase velocities lie in the velocity range between 
V, ce and K(p,a) e , which indicates that the unstable wave 
modes correspond to electron acoustic waves. 

[n] The number of spatial grid cells corresponds to N x = 
2048, and the number of velocity grid cells is chosen as N Vx = 
4096 over a velocity space that ranges from v max = 10 V, to 
v min = — 10 V, for both electrons and ions. The grid spacing 
is equal to Ax = 0.5 V tec /ui pe and the time step is equal to 
ui pe t = 0.01 . Periodic boundary conditions are imposed in the 
x direction, while open boundary conditions are imposed in 
the v A - direction. 

[ 12 ] In addition to simulations of four-component plasmas 
(Runs 0-2), we also perform a simulation with a five- 
component plasma (Run 3) consisting of a stable four- 
component plasma (Run 0) with a weak electron beam as a 
free energy source. 

3. Results 

[ 13 ] Figure 3 shows the spatial and temporal evolution of 
the electrostatic field E x component for Runs 1 and 2. The 
position is normalized by V tc Ju pe , and the time is normal- 
ized by 1 lu pe . The magnitude of electric field is normalized 

by ^ pi> V ( ( : ( jn (J c. 

[ 14 ] A wave mode with a negative phase velocity is slowly 
excited in Run 1 . The wave number at the maximum growth 
rate and the phase velocity are obtained as k x = —0.15 and 
v p ~ —0.5 V tce , respectively, from the linear analyses 
as shown in Figure 2, which is in quantitative agreement 
with the Vlasov simulation result in Figure 3 a. The excited 
wave mode seems to coalesce with the neighboring one 
and becomes isolated, as discussed in the generation of 
electron phase-space holes [e.g., Omura et al, 1996]. From 
Figure 3a, the phase velocity of the solitary structures is 


[ 10 ] In Runs 1 and 2, the temperatures of the two electron 
beams are set to be 0.5 and 0.25 of Run 0, respectively. 
The linear dispersion relations for Runs 1-2 are given in 
Figure 2. In all runs, the waves have phase velocities 
corresponding to the maximum positive gradient in the ini- 
tial velocity distribution functions (q//ov x ), which is consis- 
tent with the previous theoretical work [Singh and Lakhina, 
2001], There exists a wave with a positive growth rate only 
in the —x direction in Run 1 , while there exist waves with a 
positive growth rate in both the x and — x directions in Run 2. 
The phase velocities of unstable modes in all runs are slower 
than the thermal velocity of the core electrons V tce but faster 
than the thermal velocity of the beam electrons V t (p^ e . 
It should be noted that electron acoustic waves are generally 
defined in two-component plasmas consisting of hot and 



k x 


Figure 2. Linear dispersion relations for Runs 1-3. The 
frequency is normalized by c o pe and the wave number is 
normalized by c o p JV tec . 


3 of 7 


A05223 


UMEDA ET AL.: SIMULATION OF SOLITARY STRUCTURES 


A05223 


-4 


(a) Run 1 

-2 0 2 


xIO 

4 



(b) Run 2 



XiOpe /Vice. 


Figure 3. (a-b) Spatial and temporal evolution of electrostatic field E x for Runs 1 and 2. The position is 
normalized by V tce lu pe , and the time is nonnalized by 1 /ui pe . The magnitude of electric field is normalized 

by iOp e Vf Ce iTi e ic. 


estimated as v p ~ —0.45 V tce , which is slightly slower than 
the phase velocity of the linearly unstable mode. 

[is] The electric field profile in Figure 4a shows that the 
solitary structures have a positive electric field on the left- 
hand side and a negative electric field on the right-hand side. 
If an observer moves slower than the solitary structure, 
the polarization of the electric field E x {i) would turn from 


positive to negative, which indicates the existence of nega- 
tive potentials. It is worth showing the electron phase space 
/ ( x , v x ) and the corresponding waveforms of the electric field 
E x component, potential <fi and charge density p (Figure 4a). 
It is confirmed that there exist negative potentials and 
corresponding bipolar signatures of electric field. At the center 
of the negative potentials, the electron density is enhanced by 


(a) Run 1 



■2M 

0 , 256 512 768 1024 

X10 



XiOpe /Vtce 


(b) Run 2 



0 256 512 768 1024 



0 256 512 768 1024 

XiOpe /Vtce 


Figure 4. The electron phase space fix, v x ) and the corresponding waveforms of electric field E x , 
potential <b and charge density p at ui pe t = 4000 in (a) Run 1 and (b) Run 2. The position is normalized 
by VtcJupe, and the velocity is normalized by V tce . The magnitude of electric field is normalized by 
u pe Kce m Je, the magnitude of potential is normalized by Vj ce mje, and the magnitude of charge density 
is nonnalized by u pe mje. 
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Table 2. Characteristics of Electrostatic Solitary Structures Obtained in Runs l-3 a 


Run 


V 

w 

T 

E 

1 

0.45 V, ce ■ 

~ 1280 km/s 

150 V tce lu pe ~ 2.2 km 

150/cUp e ~ 770 /xs 

0.002 u) pe V tce m e le ~ 6.3 mV/m 

2 

0.5 V Ke - 

» 1420 km/s 

100 V lce liv pe ~ 1.5 km 

180 luj pe ~ 920 /is 

0.01 Wp 0 V t , e mJe ~ 31.5 mV/m 

3 

2.1 V tce r 

^ 5960 km/s 

50 V tce luj pe ~ 0.7 km 

2Ahjj pe ~ 122 /is 

0.0002 u ) pe V tce m e le ~ 0.6 mV/m 


a The quantities V, W, t, and E denote the velocity, pulse width, pulse duration, and electric field intensity, respectively. 


bunching of electrons in phase space. It is expected that the 
electrostatic solitary structures in Run 1 are electron acoustic 
solitary waves with negative potentials. 

[i6] In Run 2, wave modes are linearly excited in both the 
parallel and anti-parallel directions. Since waves with neg- 
ative phase velocities have higher growth rates, these waves 
show the nonlinear development into solitary structures, 
while waves with positive phase velocities show a sinusoidal 
signature. The solitary structures propagating in the anti- 
parallel direction have bipolar waveforms turning from 
negative to positive as seen in Figure 3b, which indicates the 
existence of positive potentials. The formation of solitary 
structures is due to the coalescence of potential structures 
[e.g., Omura et al, 1996]. 

[n] Figure 4b shows the electron phase space f(x, v x ) and 
the corresponding wavefonns of the electric field E x com- 
ponent and potential 6 in Run 2. It is confirmed that there 
exist electron phase-space holes with positive potentials and 
corresponding bipolar signatures of electric field. 

[is] In Figure 1, we plotted the electron velocity distri- 
bution function at the final stage (t o pe t = 4000) in Runs 1 
and 2. The final electron velocity distribution functions in 
Runs 1 and 2 show strong double peaks and are different 
from the observed velocity distribution function (Run 0). 
However, the nonlinear evolution of these two runs is very 
different. Electron acoustic solitary waves with negative 
potential exist in Run 1, while electron phase-space holes 
exist in Run 2. The final distribution function in Run 1 is 
almost the same as the initial electron velocity distribution 
function with a minor modification at v x IV tce ~ 0.5. The 
final distribution function in Run 2 shows a modification 
for v x < 0. 

[ 19 ] We can estimate the characteristics of electrostatic 
solitary structures in the Vlasov simulations with V tce ~ 
2840 km/s and u> p J2it ~ 31 kHz. Table 2 shows the char- 
acteristics of electrostatic solitary structures obtained in our 
simulation runs. We found that electrostatic solitary structures 
generated in the four-component plasma (Runs 1 and 2) 
have longer duration times (or slower propagation speeds) and 
higher wave amplitudes than the solitary structures observed 
in the magnetosheath (Event 1 in Table 3 by Lakhina et al. 
[2009]). 

[ 20 ] It should be noted that Table 3 by Lakhina et al. 
[2009] was obtained as analytic solutions to the fluid equa- 
tions for the four-component plasmas that could explain the 
characteristics of electrostatic solitary structures observed by 
the Cluster spacecraft. However, the propagation speed of 
electrostatic solitary structures obtained in their theoretical 
analysis is faster (~6000 km/s) than the drift velocities of the 
two beam components and the thermal velocities of the three 
(one core and two beam) components. Such high-speed 
electrostatic solitary structures cannot be excited in the four- 
component plasma assumed by Lakhina et al. [2009], 
because an electrostatic instability takes place at a phase 


velocity v p where the maximum positive gradient in a 
velocity distribution function is given (see the velocity dis- 
tribution functions in Figure 1 and the linear dispersion 
relations in Figure 2). Also, the electrostatic potential energy 
of solitary structures estimated by Lakhina et al. [2009] is 
much lower (1-5%) than the thermal energy of core elec- 
trons. This suggests that we might need an additional high- 
speed but weak free energy source (a positive gradient of 
velocity distribution function) at lv*l ~ 2 V tce . One of can- 
didate generation mechanisms for generating electrostatic 
solitary structures with a higher propagation speed and a 
small amplitude is the bump-on-tail instability [e.g., Omura 
et al., 1996] where electrostatic waves are excited by a high- 
speed and low-density electron beam. 

[ 21 ] In addition to Runs 1 and 2, we performed a simula- 
tion run of a five-component plasma (Run 3). The initial 
velocity distribution consists of the observed (stable) velocity 
distribution consisting of the four components (Run 0) with a 
weak electron beam as a high-speed free energy source. The 
velocity distribution becomes unstable to the well-known 
bump-on-tail instability [Omura etal., 1996]. Figure 5 shows 
the electron phase space f(x, v x ) and the corresponding 
waveforms of the electric field E x component and potential 6 
in Run 3 with the same format as Figure 4. 


V x 


E 


4 > 


Run 3 



, x 1 0 


. ■ \i i 


256 


512 


768 1024 


XUJpe / Vt ce 


Figure 5. The electron phase space fix, v x ) and the 
corresponding waveforms of electric field E x and potential 
<f> at u ) pe t = 25,000 in Run 3. The position is normalized by 
VtcJup e , and the velocity is normalized by V db . The magni- 
tude of electric field is normalized by u> pe V tce mJe, and the 
magnitude of potential is normalized by V} ce mje. 
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[ 22 ] We see that electron phase-space holes are excited in 
a narrow velocity range between —1.9 V tce and —2.2 V tce . 
The velocity distribution function is also modified in the 
narrow velocity range (see Figure 1). Such a small velocity 
space enhancement might be hardly observed by particle 
instruments [e.g., Onnirci et al., 1999]. 

[ 23 ] We estimated the characteristics of the electron 
phase-space holes (see Table 2), and found that the velocity, 
pulse width, pulse duration, and electric field of electrostatic 
solitary structures (electron phase-space holes) excited by 
the bump-on-tail instability are more similar to the analytic 
solution of the electron acoustic solitary waves obtained by 
Lakhina et al. [2009], which is also consistent with the 
amplitude and duration time of bipolar solitary pulses 
observed by the Cluster spacecraft. 

4. Summary 

[ 24 ] Using a one-dimensional electrostatic Vlasov simula- 
tion, we studied the nonlinear evolution of a four-component 
plasma consisting of one core electron component, two 
electron beams drifting in parallel and anti-parallel direc- 
tions, and one background ion component. Although the 
simulation parameters are based on Cluster observations 
made in the magnetosheath [. Lakhina et al., 2009], the para- 
meters have been slightly changed from the observations. 

[ 25 ] The present study is summarized as follows: 

[ 26 ] 1. The linear dispersion relation for the observed 
velocity distribution function (Run 0) [Lakhina et al., 2009] 
does not show positive growth rate. The Vlasov simulation 
also does not show wave growth. Thus the initial velocity 
distribution functions for the Vlasov simulation are modified 
from the observed ones in order to make them unstable. 

[ 27 ] 2. When the temperatures of the two (parallel and 
anti-parallel) electron beams are decreased slightly from the 
observed ones, the linear dispersion relation shows a posi- 
tive growth rate only in the anti-parallel direction (Run 1). 
The Vlasov simulation shows that electrostatic solitary 
structures with negative potentials are formed in the non- 
linear evolution. We expect that the solitary structures cor- 
respond to electron acoustic solitary waves with negative 
potential that propagate in the anti-parallel direction. 

[ 28 ] 3 . When the temperatures of the two electron beams 
are decreased even more, the linear dispersion relation for 
the velocity distribution function shows a positive growth 
rate in both the parallel and anti-parallel directions (Run 2). 
The Vlasov simulation shows that electrostatic solitary 
structures with positive potentials are formed in the nonlin- 
ear stage of evolution which correspond to electron phase- 
space holes [e.g., Omura et al., 1996], 

[ 29 ] 4. Duration times and wave amplitudes of electro- 
static solitary structures generated in the four-component 
plasma are longer and higher than the solitary structures 
observed by the Cluster spacecraft in the magnetosheath. 

[ 30 ] 5. Duration times and wave amplitudes of electrostatic 
solitary structures generated by the bump-on-tail instability 
are in better agreement with the Cluster observation. 

[ 31 ] It is worth noting that the characteristics of solitary 
structures generated by the bump-on-tail instability (Run 3) 
are similar to those of the electron acoustic solitary structures 
[Lakhina et al. , 2009], because phase-space holes correspond 
to phase-space structures of the weak electron beam (fifth 


component) while the other three electron components lead 
to the generation of density decreases and enhancements 
which can be described by the multifluid equations. How- 
ever, it is difficult to prove the existence of the “fifth” com- 
ponent because a small velocity space enhancement due to a 
weak electron beam might be hardly observed by particle 
instruments [e.g., Omura et al., 1999]. 

[ 32 ] Finally, it should be noted that densities and magnetic 
fields in the magnetosheath is highly turbulent [Retino et al, 
2007]. In the present study, however, we assume that the 
background ions and magnetic fields are constant because 
we focus on an electron-scale micro-instability on a time- 
scale of ca pe t ~ 25,000 which correspond to t o ci t ~ 2 in the 
magnetosheath. Large-scale fluctuations can easily change 
the types of micro-instabilities, and a large-scale multidi- 
mensional simulation is necessary. 
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